{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Classical.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OqgEsga1KFJk",
        "colab_type": "text"
      },
      "source": [
        "# Set up"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aTYBGYNxYt25",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import os, imp\n",
        "from google.colab import drive\n",
        "import numpy as np\n",
        "import h5py, os\n",
        "import cv2\n",
        "import numpy as np\n",
        "import h5py\n",
        "import gc\n",
        "from time import strftime\n",
        "from datetime import datetime, timedelta\n",
        "import matplotlib.pyplot as plt\n",
        "if not os.path.exists(\"/content/drive\"):\n",
        "    drive.mount('/gdrive', force_remount=True)\n",
        "    !ln -s '/gdrive/My Drive/GKC' drive\n",
        "\n",
        "def timestamp(info):\n",
        "    tm = datetime.now()+timedelta(hours=8)\n",
        "    print(tm.time().strftime(\"%H:%M:%S\"), info)\n",
        "os.chdir(\"/content\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "on36sEU2kxFs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def dice(pred, label):\n",
        "    return 2 * np.sum(pred * label) / (np.sum(pred) + np.sum(label))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eZhA_R5ErNnY",
        "colab_type": "text"
      },
      "source": [
        "# Threshold"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_DQsFHcAKB7M",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def gray_scale_hist(idx=1):\n",
        "    pos = []\n",
        "    fp = h5py.File(\"drive/case{}.h5\".format(idx), 'r')\n",
        "    r = np.array(fp['raw']).reshape(-1)\n",
        "    l = np.array(fp['label']).reshape(-1)\n",
        "    fp.close()\n",
        "    pos.append(r[l > 0])\n",
        "    pos = np.concatenate(pos)\n",
        "    plt.hist(pos, bins=20)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B8JV4i9fRZHU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import cv2\n",
        "from skimage.segmentation import felzenszwalb, mark_boundaries, slic, watershed\n",
        "def gray_scale_thresh(idx, thr1=77, thr2=140):\n",
        "    fp = h5py.File(\"drive/case{}.h5\".format(_), 'r')\n",
        "    r = np.array(fp['raw'])\n",
        "    l = np.array(fp['label'])\n",
        "    fp.close()\n",
        "    pred = np.zeros_like(l).astype(np.uint8)\n",
        "    a = r < thr2\n",
        "    b = r > thr2\n",
        "    pred[a * b] = 1\n",
        "    pred[:120, :, :] = 0\n",
        "    pred[:, :120, :] = 0\n",
        "    pred[:, 380:, :] = 0\n",
        "    pred[:, :, :100] = 0\n",
        "    pred[:, :, 350:] = 0\n",
        "    for i in range(120, r.shape[0]):\n",
        "        kernel = np.ones((13, 13), np.uint8)\n",
        "        pred[i, :, :] = cv2.morphologyEx(pred[i, :, :], cv2.MORPH_CLOSE, kernel, iterations=1)  # 闭运算\n",
        "        pred[i, :, :] = cv2.morphologyEx(pred[i, :, :], cv2.MORPH_OPEN, kernel, iterations=1)  # 开运算\n",
        "    return dice(pred, l)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "EkgPgmw8SYCZ",
        "colab": {}
      },
      "source": [
        "s = 0\n",
        "cnt = 0\n",
        "for _ in range(1, 11):\n",
        "    cnt += 1\n",
        "    s += gray_scale_thresh(i)\n",
        "print(s / cnt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Uv-SirkGWKMs",
        "colab_type": "text"
      },
      "source": [
        "# Classical Segmentation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Lq4kZ5ABSwMp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import skimage\n",
        "from skimage.segmentation import felzenszwalb, mark_boundaries, slic, watershed\n",
        "import os\n",
        "import matplotlib.pyplot as plt    \n",
        "import torchvision.models as models\n",
        "def T_S(idx, thre1=74, thre2=150, vis=False):\n",
        "    fp = h5py.File(\"drive/case{}.h5\".format(idx), 'r')\n",
        "    r = np.array(fp['raw'])\n",
        "    l = np.array(fp['label'])\n",
        "    fp.close()\n",
        "    pred = np.zeros_like(l)\n",
        "    for i in range(120, r.shape[0]):\n",
        "        im = r[i, 120:380, 100:350]\n",
        "        la = l[i, 120:380, 100:350]\n",
        "        pr = pred[i, 120:380, 100:350]\n",
        "        # segments = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        segments = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "        for j in range(len(segments)):\n",
        "            parser = segments == j\n",
        "            base = np.sum(parser)\n",
        "            if (base == 0):\n",
        "                continue\n",
        "            state = np.sum(im[segments == j]) / base\n",
        "            if (thre1 < state < thre2):\n",
        "                pr[parser] = 1\n",
        "    for i in range(100, 350):\n",
        "        im = r[120:, 120:380, i]\n",
        "        la = l[120:, 120:380, i]\n",
        "        pr = pred[120:, 120:380, i]\n",
        "        segments = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        #segments = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "        for j in range(len(segments)):\n",
        "            parser = segments == j\n",
        "            base = np.sum(parser)\n",
        "            if (base == 0):\n",
        "                continue\n",
        "            state = np.sum(im[segments == j]) / base\n",
        "            target = np.sum(pr[segments == j]) / base\n",
        "            if (74 < state < 150 and target > 0.9):\n",
        "                pr[parser] = 1\n",
        "            elif target < 0.5:\n",
        "                pr[parser] = 0\n",
        "    if (vis):\n",
        "    \n",
        "        fig = plt.figure()\n",
        " \n",
        "        ax1 = fig.add_subplot(231)\n",
        "        ax1.imshow(r[250, :, :])\n",
        " \n",
        "        ax1 = fig.add_subplot(232)\n",
        "        ax1.imshow(l[250, :, :])\n",
        " \n",
        "        ax1 = fig.add_subplot(233)\n",
        "        ax1.imshow(pred[250, :, :])\n",
        " \n",
        "        ax1 = fig.add_subplot(234)\n",
        "        ax1.imshow(r[:, :, 250])\n",
        " \n",
        "        ax1 = fig.add_subplot(235)\n",
        "        ax1.imshow(l[:, :, 250])\n",
        " \n",
        "        ax1 = fig.add_subplot(236)\n",
        "        ax1.imshow(pred[:, :, 250])\n",
        "    return dice(pred, l)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JqxvPTqj-PC1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "s = 0\n",
        "cnt = 0\n",
        "for _ in range(1, 11):\n",
        "    s += T_S(_)\n",
        "    cnt += 1\n",
        "    print(s)\n",
        "print(s / _)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ze153G6eUae4",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import skimage\n",
        "from skimage.segmentation import felzenszwalb, mark_boundaries, slic, watershed\n",
        "import os\n",
        "import matplotlib.pyplot as plt    \n",
        "import torchvision.models as models\n",
        "def S_T(idx, thre1=74, thre2=150, vis=False):\n",
        "    fp = h5py.File(\"drive/case{}.h5\".format(idx), 'r')\n",
        "    r = np.array(fp['raw'])\n",
        "    l = np.array(fp['label'])\n",
        "    fp.close()\n",
        "    \n",
        "    pred = np.zeros_like(l)\n",
        "    for i in range(100, 350):\n",
        "        im = r[120:, 120:380, i]\n",
        "        la = l[120:, 120:380, i]\n",
        "        pr = pred[120:, 120:380, i]\n",
        "        segments = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        # segments = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "        for j in range(len(segments)):\n",
        "            parser = segments == j\n",
        "            base = np.sum(parser)\n",
        "            if (base == 0):\n",
        "                continue\n",
        "            state = np.sum(im[segments == j]) / base\n",
        "            if (thre1 < state < thre2):\n",
        "                pr[parser] = 1\n",
        "    for i in range(120, r.shape[0]):\n",
        "        im = r[i, 120:380, 100:350]\n",
        "        la = l[i, 120:380, 100:350]\n",
        "        pr = pred[i, 120:380, 100:350]\n",
        "        segments = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "        # segments = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        for j in range(len(segments)):\n",
        "            parser = segments == j\n",
        "            base = np.sum(parser)\n",
        "            if (base == 0):\n",
        "                continue\n",
        "            state = np.sum(im[segments == j]) / base\n",
        "            target = np.sum(pr[segments == j]) / base\n",
        "            if (74 < state < 150 and target > 0.9):\n",
        "                pr[parser] = 1\n",
        "            elif target < 0.5:\n",
        "                pr[parser] = 0\n",
        "    \n",
        "    if (_ == 2):\n",
        "        \n",
        "        fig = plt.figure()\n",
        "        im = r[250, :, :]\n",
        "        segmentsF = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        segmentsS = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "\n",
        "        ax1 = fig.add_subplot(221)\n",
        "        ax1.imshow(mark_boundaries(im, segmentsF))\n",
        " \n",
        "        ax1 = fig.add_subplot(222)\n",
        "        ax1.imshow(mark_boundaries(im, segmentsS))\n",
        "        \n",
        "        \n",
        "        im = r[:, :, 250]\n",
        "        segmentsF = felzenszwalb(im, scale=5000, sigma=0.01, min_size=500)\n",
        "        segmentsS = slic(im, n_segments=150, compactness=0.1, sigma=1)\n",
        "        ax1 = fig.add_subplot(223)\n",
        "        ax1.imshow(mark_boundaries(im, segmentsF))\n",
        " \n",
        "        ax1 = fig.add_subplot(224)\n",
        "        ax1.imshow(mark_boundaries(im, segmentsS))\n",
        "    return dice(pred, l)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2peG0qed2OHl",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "s = 0\n",
        "cnt = 0\n",
        "for _ in range(1, 11):\n",
        "    s += dice(pred, l)\n",
        "    cnt += 1\n",
        "print(s / cnt)"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}